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COMBINED REFLECTION AND DIFFRACTION BY A VERTICAL WEDGE 


PART I: INTRODUCTION 


1. The boundary value problem of linear wave reflection and diffraction 
by a vertical wedge of arbitrary wedge angle has been well formulated and 
presented by Stoker (1957) among many other investigators. The technique to 
obtain an analytical solution for the problem is also depicted in the cited 
book. However, analytical solutions are not available for the problem, except 
for the special case of wave diffraction by a thin semi-infinite breakwater, 
that is, a wedge with wedge angle equal to zero. 

2. The solution of the thin semi-infinite breakwater was presented in 
the dimensionless diffraction diagrams by Wiegel (1962). The diagrams have 
been especially useful in preliminary engineering design and have been in- 
cluded in the Shore Protection Manual (SPM) (1984). Although equally useful, 
the combined reflection and diffraction diagrams are not available, perhaps 
because of the complexity of the diagrams which makes them difficult to create 
without using modern high-speed computers for computation and graphing. 

3. The objectives of the present study are (a) to obtain an analytical 
solution for the combined wave reflection and diffraction by a vertical wedge 
of arbitrary wedge angle subject to excitation of a plane simple harmonic wave 
train coming from infinity and (b) to provide the combined reflection and dif- 
fraction diagrams. The diagrams included in this report have two cases: one 
for a thin semi-infinite breakwater and the other for a 90-deg vertical wedge. 
Subroutine WEDGE for computing the combined reflection and diffraction by a 
vertical wedge of arbitrary wedge angle is also documented in the report 
(Appendix A). 


PART II: BOUNDARY VALUE PROBLEM 
Mathematical Formulation 


4. In this study our primary interest is the wave reflection and 
diffraction by a vertical wedge of arbitrary wedge angle in a constant water 
depth h * subject to the excitation of monochromatic incident waves of infi- 
nitesimal amplitude coming from infinity. Let (r,0,z) be cylindrical coordi- 
nates, with z= 0 representing the undisturbed water free surface and upward 
direction representing the positive z-axis. The tip of the wedge is chosen to 
be the origin of the coordinates and two rigid walls of the wedge to coincide 
with 6=0 and 6= 85 » respectively, as illustrated in Figure 1. Carte- 
sian coordinates (x,y,z), corresponding to the cylindrical coordinates, are 


also occasionally used and shown in the same figure. Therefore, the wedge 


zZ 


Figure 1. A vertical wedge of 
arbitrary wedge angle 


* For convenience, symbols and abbreviations are listed in the notation 
(Appendix B). 


angle is 27 - on » and the water region is defined by oF 2620 and 
0222-h. 

5. The velocity field for the wave reflection and diffraction in an 
ideal fluid can be represented by the velocity potential function 4(r,®,z,t) 
which must satisfy the Laplace equation, where t is the temporal coordinate. 
We assume that the waves are sinusoidal in time with radian frequency w . 
Water depth is constant, and the bottom is rigid and impermeable. Therefore, 
the vertical and temporal components of the velocity potential function, which 
follow from separation of variables, can be factored out and the velocity po- 


tential written as 


mn cosh k(z + h) iwt 
O(r,8,z,t) = Ay cosh kh o(r,8)e (1) 
where 
Ay = -iga,/w 
i= yal 
g = gravitational acceleration 
aS incident wave amplitude 


k = wave number 

¢ = horizontal component of the velocity potential function 

6. Substituting Equation 1 into the Laplace equation and using both 
the kinematic and dynamic boundary conditions at the free surface, the Laplace 
equation is then reduced to the Helmholtz equation which is written in polar 


coordinates as follows: 


ES aot gine ai gh kro =0 (2) 


where k must be a real number and satisfy the dispersion relationship 
2 
W = gk tanh kh (3) 


7. The free surface displacement NN from the mean water level z= 0 
can be obtained from linear wave theory and is represented as 
1 ao 


= Oe iwt 
n(r,8,t) = a 8G aj o(r, Oe (4) 


8. Thus only the horizontal part of the velocity potential function 9 
is needed to be determined as a solution of Equation 2 in the water region 
on > 620 , with the following boundary conditions at the rigid and imperme- 


able walls of the wedge: 


ad _ 0 


ona ae @ = 0 ernal 85 (5) 


9. A condition at infinity is also required to ensure a unique solu- 
tion. The classic approach is to use the Sommerfeld radiation condition at 
infinity which states that the scattered wave 0 must behave like a cylin- 


drical outgoing progressing wave at infinity such that ' 


a 
lim vr (ste + sn4,) = 0 (6) 


The total wave represented by 9 is the linear superposition of an incident 
wave %; » a reflected wave from the the ®=0 wall of the wedge o >» and 


the scattered wave Os from the tip of the wedge. 


9= 6,4 6,4 9, (7) 
Equation 6 can be satisfied if 
eo tkr 
od ~ at ro (8) 
= ise 


10. The incident wave coming from a large distance from the tip of the 
wedge is assumed to be a plane progressive wave of amplitude ay and incident 
angle a to the x-axis as given by 


=f etkr cos (6-a) 


(9) 


Consequently, the perfectly reflected wave from the y = 0 wall of the wedge 


is 


a etkr cos (6+0) (10) 


Thus the boundary value problem (in which the governing equation is Equa- 
tion 2, the boundary condition is Equation 5, and the radiation condition is 


Equation 6) is completely formulated. 


Analytical Solution 


11. Analytical solution to the problem formulated in the preceding 
section is obtained by following the solution technique by Stoker (1957). To 
obtain the solution, the water region is divided into three subregions--I, II, 
and III--by the incident wave ray passing through the tip of the wedge and the 
reflected wave ray reflected away from the tip of the wedge, as shown in Fig- 
ure 2. Obviously, the total wave in subregion I is the sum of the incident, 
reflected, and scattered waves; the total wave in subregion II, where the 
reflected wave does not exist, is the sum of the incident and scattered waves; 
and the total wave in subregion III, where the incident and reflected waves 


have been shaded out, is only the scattered wave. For certain combinations of 


y 


Tt 


Figure 2. Three subregions and the wedge 


the wedge angle and incident wave angle, subregions II and III may not exist 


at all. In general, the solution function can be written as 


$= ¢(r,8) + 6 (r,8) (11) 
where 
5 + om T-a>O>0 
o (F298) = o; Tta>O6@>7T-a (12) 
0 oe > @ > ta © 


The equation reveals that oF is the sum of the incident and reflected waves 
5 and $. and is a known function. The scattered wave oe is the only un- 
known function to be determined in the problem. Nevertheless, the total wave 
¢ instead of the scattered wave 0. is the desired solution to be obtained 
in this study. 

12. The solution for the total wave $ is pursued. The finite cosine 


transform of 9 , denoted by $ » is introduced by the formula 


VT 


o(kr,n) = i o(kr,®) cos ut dé (13) 


where n= 0,1,2,... are integers, and v is related to the wedge angle as 


defined by 


@ = vt (14) 


Applying the finite cosine transform and using the boundary condition in Equa- 


tion 5, Equation 2 becomes 


etree +r = + (IED - (2) ¢=0 (15) 


Equation 15 is a form of the Bessel equation for which general solutions are 
Ja sy 6k) and Y/vike) ; 


respectively. Since Ya/vik®) are singular at the origin, the solution is 


the Bessel functions of the first and second kinds, 


chosen to be 


o(kr,n) = a Jy 6k) (16) 


where a, are constants to be determined. 
13. Taking the finite cosine transform of Equation 11 and using Equa- 


tion 16, we have 


VT VOT 


a Jy (kX) - i 6 cos = dé (17) 
0 0 


~o- 
n 
Q 
ie} 
n 
Is 
Q 
@D 
i] 


or 


(kr) - co (18) 


2 | 


Son ss 


Then applying the operation lim... vr(9/dr + ik) to both sides of Equa- 


tion 18, and using the Sommerfeld radiation condition (Equation 6) we have 


VT 
F 0 , ng “ 
Lim (Fe is ix) a. Tes) © i 6, cos * ae| = 0 (19) 
( 


14. Equation 19 can be asymptotically evaluated to determine ale 
Firstly, the first term involving the Bessel function is evaluated. The func- 
tion J py Sk) at reo behaves asymptotically (Abramowitz and Stegun 1964) 


as follows: 


Jay (ke) ~ Vo cos (ex - = - z) (20) 


Hence, we have 


Lim "(Fe a, th) py : E: ei (kr-nn/2vtn/4) (21) 


Secondly, the second term involving the integral of oF is evaluated. The 
asymptotic behavior of the integral over 9% = (0,v™) and at large distance 
r>© can be found by the method of stationary phase. The integral, after sub- 


stituting oP from Equations 9, 10, and 12, can be written as 


vt TO. 
e ie ° 6 
i d cos ng dé -f E: Sees=e)) + epliee cos(842) | Coste 
fo) v av) 
0 0 
TO. 
i C= () 
+ gue eerie) cos = do (22) 
T-a 


In the integrals, there are three points of stationary phase at 9 = 4% and 
6=7 +a. If the same argument as that of Stoker (1957) is followed, of the 
three contributions only the first one 9 =a furnishes a nonvanishing con- 
tribution for r*@ when the operator Y“r(9/dr + ik) is applied to it. The 
physical significance of this statement is that only the incident wave is 


effective in determining the cosine coefficients of the solution. Therefore, 


VT 
lin “e(% + ax) f 9, cos ne d6 ~ 2V2Tk cos = sare) (23) 
0 


Substituting Equations 21 and 23 into Equation 19, we obtain the unknown co- 


efficients a 


na eint/2v 


anes ZECOS Ra (24) 
15. Since the solution ¢ in the cosine series expression is 
ee ee 2 = nO 
6(r,8) = = B(r,0) += JY F(xsn) cos & (25) 
n=1 


10 


the solution is obtained by substituting Equations 16 and 24 into Equation 25 


as follows: 


inn/2v a. ng 


o(r,9) = = JQ (kr) + 2 » e Jay ke) cos - cos —— (26) 


n=1 


Equation 26 is the solution for the combined wave reflection and diffraction 
by a vertical wedge of arbitrary wedge angle and is considered to be extended 
from the solution by Stoker (1957) who only solved the problem of a thin semi- 
infinite breakwater. The solution in Equation 26 and the one by Stoker are 
not only in nonclosed form but also in terms of Bessel functions. It seems 
that the calculations of the solutions are very difficult without using a mod- 
ern high-speed computer. This is probably the reason why Stoker arrived at 
his solution expressed in the same cosine series but did not use it to calcu- 
late the result. Instead, he further transformed the expression into a very 
complex integral form for further approximation in calculating the result. 

16. Notably, the solution at the origin point is obtained by simply 


substituting r= 0 into Equation 26 to arrive at 
2 
$(0,8) == (27) 


Therefore, wave response at the origin point depends only on the wedge angle 


and does not depend on the incident wave angle. 


Two Special Cases 


17. The solutions for two special cases are used to verify Equation 26: 
one for the case of a thin semi-infinite breakwater and the other for the case 
of an infinite wall extending from x = -~ tom, 

18. The vertical wedge should reduce to a thin semi-infinite breakwater 
as the wedge angle reduces to 0 deg. Therefore, solution of the combined wave 
reflection and diffraction by a thin semi-infinite breakwater is obtained by 
substituting v=2 (that is, oR = 27) into Equation 26 which then becomes 


11 


OGe,@)) = JQ (kr) + 2 » einm/4 5 (kr) cos = cos x (28) 


n=l 


n/2 


Equation 28 is precisely the same one obtained by Stoker (1957). 

19. The vertical wedge should also become an infinite wall extending 
from x = -© to © with the water occupying only the half plane of y20 as 
the wedge angle increases to 180 deg. In this situation the scattered wave is 
absent from the solution, and the total wave is only the sum of the incident 


and reflected waves as follows: 


GeO) = eikr cos (6-0) + eikr cos (8+0) ag 


After expansion of the exponential functions in terms of Bessel functions 


(Abramowitz and Stegun 1964), Equation 29 becomes 


foe) 


$(r,8) = 2 Jo (kr) ap 2 » i"J (kr) cos n@ cos n@ (30) 


n=1 


Equation 30 is the same equation reduced from Equation 26 by substituting 


v= 1 into it. 


12 


PART III: CALCULATION AND RESULTS 


20. Results of the combined reflection and diffraction by a wedge of 
arbitrary wedge angle can be calculated from Equation 26. Since the solution 
is not only in terms of Bessel functions but also in a nonclosed form, the 
computer program WEDGE is therefore written to calculate the solution. 

21. In the program the subroutine BESJ for calculating Bessel function 
of fractional or integer order was used. The subroutine was originally writ- 
ten by Amos, Daniel, and Weston in 1975 (Morris 1984) and is collected in the 
Naval Surface Weapons Center Library of Mathematics Subroutines (Morris 1984). 

22. In the calculation the summation of the infinite terms in Equa- 
tion 26 was carried out to the term which is preceded by eight successive 
terms of the absolute value of the Bessel function, all equal to or less than 
Oman The solution has a truncation error less than Omen and it is of the 
order of one. 

23. In this study, results of the combined wave reflection and diffrac- 
tion for the wedge are calculated for two cases: one for a vertical wedge of 


O-deg wedge angle and the other for a vertical wedge of 90-deg wedge angle. 


Vertical Wedge of 0-Deg Wedge Angle 


24. When the wedge angle is equal to zero, the wedge is actually a thin 
semi-infinite breakwater extending from x =0 to ». Figure 3 shows the 
thin semi-infinite breakwater along with the polar coordinates. In this case 
the diffraction results for various incident wave angles in the water region 
from 6 =m to 2m and r/dA < 10 , where A is the incident wave length, 
have already been presented by Wiegel (1962) and are shown in the SPM, Vol- 
ume I (1984). The present results combine reflection and diffraction effects 
and cover the water region from 6 =0 to 2m and r/dA < 10. Therefore, 
the present results for this particular case can be considered to be a comple- 
mentary and extended version to the ones in the SPM. 

25. In this study wave response was calculated at 1,460 grid points 
intersected by r/\ = 0.5,(0.5),10 , which means that the values of r/) are 
from 0.5 to 10.0 with each value increment being 0.5. Hereafter, all similar 


expressions are to be interpreted in the same way (e.g., 96 = 0,(7/36),2m for 
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Figure 3. Thin semi-infinite breakwater and polar coordinates 


the incident wave angle a = 0,(mn/12),7) . The wave response at the origin 


point is obtained substituting v = 2 into Equation 27, as follows: 


o(0,6) = 1 (31) 


Those calculated values were used to interpret the value for each non- 


overlapping pixel of size O.lr/A by O.lr/X in the area within the 10r/A 


radius from the origin. A diagram was then constructed by patching those 


The wave response diagrams for each incident 


pixels over the entire area. 
Notably, the values in the dia- 


wave angle are shown in Figures 4 through 15. 


grams constitute the amplification factor which is defined as the ratio of the 


total wave height to the incident wave height. Therefore, in subregions II 
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Figure 4. Amplification factor diagram for the thin semi-infinite 
breakwater for incident wave angle = 0 deg 


and III (as defined in Figure 2) where the reflected wave is absent, the am- 
plification factor is essentially the diffraction coefficient as defined in 
the SPM. 

26. Figures 4 through 15 reveal that the amplification factors in sub- 
region I change very rapidly between 0 and 2.35 over the subregion, and the 
diagram patterns become very complex because of the interesting superposition 
of the incident, reflected, and scattered waves. (In the legend of Figures 4 
through 15, the width of the pixel is one incident wave length, and the values 
are amplification factors.) Such patterns would be very difficult to con- 
struct without using a high-speed computer and computer graphics. In sub- 
regions II and III, the amplification factors change smoothly from 1.15 


roughly along the reflected wave ray reflected from the origin point to nearly 
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Figure 7. 


Amplification factor diagram for the thin semi-infinite 
breakwater for incident wave angle = 45 deg 
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Figure 13. Amplification factor diagram for the thin semi-infinite 
breakwater for incident wave angle = 135 deg 
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0.00 at the back wall of the wedge in the shadow zone. The diagram patterns 
are relatively smooth and simple. 

27. Notably, the diagrams do not include phase information of the wave 
response which is usually unimportant in most engineering practice. Should 
the phase of the wave response need to be known, one can always use the com— 
puter program WEDGE to calculate it. 

28. The contour diagrams of the amplification factor, similar to the 
ones presented by Wiegel (1962) and shown in the SPM (1984), were also plotted 
as typically shown in Figure 16. Examination of those contour diagrams indi- 
cates that, in subregion III, the results are identical to Wiegel's results. 
But in subregion II the contour patterns for the amplification factor (or the 
diffraction coefficient K" used in the SPM (1984)) equal to 1.0; thus the 
present results are far more complicated than Wiegel's. The author believes 
that Wiegel's results may lose accuracy because of insufficient resolution of 
the computational tools during the late fifties and early sixties. Neverthe- 
less, such inaccuracies are usually either tolerable or immaterial in most 
engineering practice. 

29. The contour patterns in subregion I are very complex, and it is 
difficult to track specific contours. Therefore, for clarity only the patched 


diagrams are presented, and the contour diagrams are omitted in this report. 
Vertical Wedge of 90-Deg Wedge Angle 


30. When the wedge angle is equal to 7/2(6, = 31/2) , the vertical 
wedge occupies the entire fourth quadrant of the space as shown in Figure 17. 
Wave response was calculated at 1,100 grid points intersected at r/d 
= 0.5,(0.5),10.0 and 6 = 0,(1/36),31/2 for the incident wave angle a 
= 0,(1/12),7 . The wave response at the origin is obtained by substituting 


v = 1.5 into Equation 27 as follows: 


6(0,0) =< (32) 


Those calculated results were used to construct the amplification factor dia- 
grams by following the same procedures described for the case of the thin 
semi-infinite breakwater. The diagrams are shown in Figures 18 through 27. 
Because of symmetry of the results, the diagrams for the incident wave angle 


a > 31/4 can be obtained from those for an incident wave angle of tT-a. 


2H, 
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Figure 16. Amplification factor contour diagram for the thin 


semi-infinite breakwater for incident wave angle = 90 deg 
Therefore, the diagrams for the incident wave angles 
omitted in this report. 


Sil, 


a = 1/6,(1/12),7 are 


The diagrams indicate that, for each corresponding incident wave 


angle, the results in subregion I are very similar to those obtained for the 


vertical wedge of O0-deg wedge angle. However, the results in subregions II 


and III from both wedges are discernibly different. 
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A 90-deg wedge and polar coordinates 
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Figure 17. 


Figure 18. Amplification factor diagram for the 


90-deg wedge for incident wave angle 
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Figure 19, Amplification factor diagram 
90-deg wedge for incident wave angle = 
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Figure 21. Amplification factor diagram for the 
90-deg wedge for incident wave angle = 45 deg 
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Figure 22. Amplification factor diagram for the 


90-deg wedge for incident wave angle 
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Figure 24. Amplification factor diagram for the 
90-deg wedge for incident wave angle = 90 deg 
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PART IV: CONCLUSION 


32. An analytical solution for the combined wave reflection and dif- 
fraction by a vertical wedge of arbitrary wedge angle is obtained and ex- 
pressed in Equation 26. The analytical solution is in terms of Bessel 
functions and in nonclosed form. The computer subroutine WEDGE, written for 
calculating the solution, is documented in Appendix A. 

33. The amplification factor diagrams for a vertical wedge of O-deg 
wedge angle and a vertical wedge of 90-deg wedge angle are calculated and pre- 
sented. The calculated results indicate that the wave response in subre- 
gion I, where the incident, reflected, and scattered waves all exist, is in a 
very complex pattern with the amplification factor varying from 2.35 to 0.0 
over the subregion. The wave response in subregions II and III is a rela- 
tively simple pattern with the amplification factor decreasing from 1.15 
roughly along the reflected wave ray reflected from the origin point to nearly 
0.00 at the back wall of the wedge in the shadow zone. 

34. Diagrams of the special case of a vertical wedge of O-deg wedge 
angle can be considered complementary and extended versions to the ones pre- 


sented in the SPM (1984). 
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APPENDIX A: SUBROUTINE WEDGE 


1. Subroutine WEDGE is used to calculate the value of 94 in Equa- 

tion 26 which is generally a complex number. Its absolute value is the ampli- 
fication factor, and its phase is the phase indicator from the phase of the 
incident wave. As mentioned in the main text, 4 is a function of the Bessel 
function of either fractional or integer order, depending on the wedge angle, 
and is the aummaeien of a series of infinite terms. The subroutine BESJ, doc- 
umented in the Naval Surface Weapons Center (NSWC) Library of Mathematics Sub- 
routines (Morris 1984)* is used in the WEDGE subroutine. The programming of 
the WEDGE subroutine is very straightforward if a truncation term in the se- 
ries in Equation 26 is determined. The program is written in FORTRAN language 


and is listed in this appendix. 


Description 


2. The following subroutine is available for computing 9 in 


Equation 26: 
CALL WEDGE(F, FABS, FPHA,XRL, XTH, WEDGEA, WAVEA, IDX) 


where the arguments are all real values except F which is a complex value. 


Input arguments are as follows: 


a. (XRL,XTH)=(r/A,60) where (r,9) are polar coordinates of the lo- 
cation where $ is to be computed, and A is the incident wave 
length. Therefore, XRL is the radius vector or radius distance 
normalized by the incident wave length. XTH is the vectorial 
angle in degree. 


- WEDGEA = wedge angle in degree. 


Io 


c. WAVEA = incident wave angle in degree. 
d. IDX = an index (set to 0 in this subroutine). 


Output arguments are as follows: 


* References cited in the Appendix can be found in the References at the end 
of the main text. 


Al 


a. F= @ in Equation 26 (wave response normalized by incident 
wave amplitude). 


b. FABS = amplification factor, the absolute value of 9. 


ig 
ies] 
Las] 
Il 


phase difference, the phase of $6. 


Example and Test Run 


3. To serve as an example as well as to ensure the subroutine is the 
correct one, the user should run the test program listed in Figure Al and make 


sure the output is the same as that listed in Table Al. 


FROGRAM WEDGE! (INPUT ,OUTPUT, TAPES=INPUT, TAFES=QUTPUT) TESTI 
CONPLEX F TEST2 
PI=3. 141592654 TESTS 
WRITE (4,4) TEST4 

4 FORMAT(//3%,' WEDG-ANG WAV-ANG LOCATION WAVE RESPONSE TESTS 
{ ABS-VAL PHASE’/3X,’ (DEG) (DEB) XRL XTH(DEB) (NORMTESTS 
PALIZED) AMP-FAC (RAD) ‘/3X,35¢° -*)/) TEST? 

WEDGEA=90, TEST 
WAVEA=135. TESTS 
XRL=2.0 TESTIO 
¥TH=30. TEST11 
IDY=0 TESTI2 
CALL WEDGE (F,FABS,FPHA, XRL,XTH,WEDGEA,WAVEA, IDX) TESTI 
WRITE(6,40) WEDGEA,WAVEA,XRL,XTH,F,FARS,FPHA TESTIA 
40 FORMAT(1X,8F9. 2) TESTIS 
STOP TESTI4 
END TESTI7 


Figure Al. Computer program list 1 


Table Al 


Sample Output of the Test Program 


WEDG-ANG WAV-ANG LOCATION WAVE RESPONSE ABS-VAL FHASE 
(DEG) (DEG) RRL XTH(DEG) (NGRMALIZED) AMP-FAC iRAD) 
90.00 135.00 2.00 40.66 a5 Ail 268 a fe Bo Mil 


A2 


4. The input arguments in the test programs are as follows: 
a. WEDGEA = 90; the wedge angle is 90 deg. 
b. WAVEA = 135; the incident wave angle is 135 deg. 


c. XRL = 2.0 and XTH = 30; the location of the wave response $4 to 
be computed is at radial vector of two incident wave length dis-— 
tances and vectorial angle of 30 deg in polar coordinates. 


d. IDX = 0; the index IDX is set to 0. 
This case is shown in Figure A2. The outputs are given in Table Al which is 
self-explanatory. 

5. If the location of the wave response $ to be computed is a very 
large distance, for example XRL greater than 18 or so, the output might print 
the message that the number of terms is insufficient for summation in comput- 
ing $ . In this situation, the user must replace the integer 200 in 
PARAMETER(NN=200), listed in Card WEDGE15 in the SUBROUTINE WEDGE list, by a 


larger integer to ensure the accuracy. 


Figure A2. Example problem 
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Subroutine WEDGE Listing 


6. Subroutine WEDGE, which is listed in this section, calls subroutine 
BESJ which, in turn, calls subroutine JAIRY and function GAMLN. Subroutines 
BESJ and JAIRY and function GAMLN are borrowed from the NSWC Library of Mathe- 
matics Subroutines (Morris 1984). Including these borrowed subroutines here 
in the list is only for the purpose of allowing subroutine WEDGE to be self- 
contained and complete. The computer program of subroutine WEDGE is listed in 


Figure A3. 
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SUBROUTINE WEDGE(F,FABS,FPHA,XRL,XTH,WEDGEA,WAVEA, IDX) 


THIS COMPUTER PROGRAM WAS WRITTEN BY H.S. CHEN OF CERC IN 
1986 UNDER CORP’ CIVIL WORK R&D PROGRAM. NEITHER ANY OF 
AGENCIES NOR ANY INDIVIDUAL ASSUMES ANY LEGAL LIABILITY OR 
RESPONSIBILITY FOR THE ACCURACY OF THE PROGRAM. 


Es I 2 


WAVE REFLECTION AND DIFFRACTION BY A VERTICAL WEDGE. 


INPUT : XRL = R/L, (RADIUS VECTOR) / (WAVE LENGTH) 
XTH = VECTORIAL ANGLE IN DEGREE 
WEDGEA = WEDGE ANGLE IN DEGREE 
WAVEA = INCIDENT WAVE ANGLE IN DEGREE 
OUTPUT : F = FHI, WAVE RESPONSE NORMALIZED BY INCIDENT 
WAVE AMPLITUDE 
FABS = ABSOLUTE VALUE OF PHI, THE AMPLIFICATION 
FACTOR 
FPHA = PHASE OF PHI IN RADIAN. 


PARAMETER (NN=200) 

DIMENSION BJ (NN) ,W(NN) , XM(NN) 

COMPLEX F,1M 

DATA TOLR/1.E-8/,1TER/8/ 

PI=3.141592654 

CPI=PI/180. 

XKR=XRL#2. O#P I 

TH=XTH#CPI 

WA=WAVEA*CPI 

XNU=(360.-WEDGEA) /180. 

XM(1)=1.0 

IF(IDX.NE.0) GOTO 14 

CALL BESJ(XKR,0.0,1,W,NZ) 

Bd(1)=W(1) 

ICOUNT=0 

DO 10 N=1,NN 

NL=N+1 

ICOUNT=ICOUNT+1 

IF(ICOUNT.LE.ITER) GOTO 8 

NNN=N 

GOTO 14 
8 XM(N1)=FLOAT(N) /XNU 

M=INT(XM(NL) ) 

ALPHA=XM(NI) -M 

ML=M+1 

CALL BESJ(XKR,ALPHA,M1,W,NZ) 

IF(N1.EQ.NN) WRITE(4,9) XKR,ALPHA,ML,W(ML) ,NZ 
9 FORMAT(/'’ ##** NO. OF TERMS FOR SUMMATION IS INSUFFICIENT *#*%', 

1 /’ XKR,ALPHA,M1,W(ML),NZ =',2F10.4,15,E15.6,15/) 

BJ (N41) =W(ML) 


Figure A3. Computer program list 2 (Sheet 1 of 25) 
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WEGDEL 


Re HE EF € RE KR HE KE KK KR RR KR K€ EK KK HR HE KR KE KE RE & KF KF & E¥NOTICEL] 


*¥NOTICE2 
*NOTICES 
*NOTICE4 
*¥NOTICES 


%¥ # FF HK He EH EK HK KF RH HE KH HF HH HH HK KR HR HK HK HH HE HE KR KF RNOTICES 


-WEDGE2 
WEDGES 
WEDGE4 
WEDGES 
WEDGES 
WEDGE7 
WEDGES 
WEDGE? 
WEDGE10 
WEDGE1L1 
WEDGE12 
WEDGELS 

-WEDGE14 
WEDGE1S 
WEDGE14 
WEDGE17 
WEDGE18 
WEDGE19 
WEDGE20 
WEDGE21 
WEDGE22 
WEDGE23 
WEDGE24 
WEDGE25 
WEDGEZ4 
WEDGE27 
WEDGE28 
WEDGE29 
WEDGES 
WEDGES1 
WEDGES2 
WEDGE3S3 
WEDGES4 
WEDGESS 
WEDGE364 
WEDGES7 
WEDGESS8 
WEDGE39 
WEDGE40 
WEDGE41 
WEDGE42 
WEDGE43 
WEDGE44 


oa 
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IF (ABS(BJ(N1)).GT.TOLR) ICOUNT=0 
10 CONTINUE 
14 CONTINUE 
F=BJ(1)/2. 
DO 20 N=1,NNN 
N1=N+1 
XMN=XM (NL) 
TM= (0.0, 1.0) ##XMN#BJ (NI) #COS (XMN#WA) #COS (XMN*TH) 
F=F+TM 
20 CONTINUE 
F=4./XNU*F 
FR=REAL (F) 
FI=AIMAG(F) 
FABS=SORT (FR¥FRtFI¥#FI) 
IF (WA.LE.1.E-8) FABS=FABS/2. 
IF(FABS.LT.TOLR) GOTO 30 
FPHA=ATAN2(FI,FR) 
RETURN 
30 FPHA=0.0 
RETURN 
END 
SUBROUTINE BESJ(X,ALPHA,N,Y,NZ) 
WRITTEN BY D.E. AMOS, S.L. DANIEL AND M.K. WESTON, JANUARY, 1975. 
REFERENCE SAND-75-0147 


ABSTRACT 
BESJ COMPUTES AN N MEMBER SEQUENCE OF J BESSEL FUNCTIONS 
J/SUB(ALPHAtK-1)/(X), K=1,...,N FOR NON-NEGATIVE ALPHA AND Xx. 
A COMBINATION OF THE POWER SERIES, THE ASYMPTOTIC EXPANSION 
FOR X TO INFINITY AND THE UNIFORM ASYMPTOTIC EXPANSION FOR 
NU TG INFINITY ARE APPLIED OVER SUBDIVISIONS OF THE (NU,X) 
PLANE. FOR VALUES OF (NU,X) NOT COVERED BY ONE GF THESE 
FORMULAE, THE GRDER IS INCREMENTED OR DECREMENTED BY INTEGER 
VALUES INTO A REGION WHERE ONE OF THE FORMULAE AFPLY. BACKWARD 
RECURSION IS APPLIED TO REDUCE ORDERS BY INTEGER VALUES EXCEPT 
WHERE THE ENTIRE SEQUENCE LIES IN THE OSCILLATORY REGION. IN 
THIS CASE FORWARD RECURSION IS STABLE AND VALUES FROM THE 
ASYMPTOTIC EXPANSION FOR X TO INFINITY START THE RECURSION 
WHEN IT IS EFFICIENT TO DO SO. LEADING TERMS OF THE SERIES AND 
UNIFORM EXPANSION ARE TESTED FOR UNDERFLOW. IF A SEQUENCE IS 
REQUESTED AND THE LAST MEMBER WOULD UNDERFLOW, THE RESULT IS 
SET TO ZERO AND THE NEXT LOWER ORDER TRIED, ETC., UNTIL A 


MEMBER COMES ON SCALE OR ALL MEMRERS ARE SET TO ZERO. OVERFLOW 
CANNOT OCCUR. BESJ1 CALLS SUBROUTINE JAIRY AND FUNCTION GAHMLN. 
DESCRIPTION OF ARGUMENTS 
INPUT 
X | Hollen @ 


ALPHA - ORDER OF FIRST MEMBER OF THE SEQUENCE, ALPHA.GE.90 


N - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1 
OUTPUT 
Y - A VECTOR WHOSE FIRST N COMPONENTS CONTAIN 


VALUES FOR J/SUB(ALPHAtK-1)/(X), K=1,...,N 


Figure A3. (Sheet 2 of 25) 
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BESJ114 
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BESJ116 
BESJ117 
BESJ118 
BESJ119 
BESJ120 
BESJ1i21 
BESJ1i22 
BESJ1253 
BESJ124 
BESJ125 
BESJ126 
BESJ127 
BESJ128 
BESJ129 
BESJ130 
BESJi31 
BESJ132 
BESJ133 
BESJ134 
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NZ - ERROR 
NZ=0 
NZ=-1 
NZ=-2 
NZ=-3 
NZ.GT. 


ERROR CONDITIONS 


IMPROPER INPUT AR 
UNDERFLOW 


DOUBLE PRECISION DX,T 
DIMENSION Y(N) 

DIMENSION Ciii,10) 
DIMENSION C1(88) 
DIMENSION A1(952), 
DIMENSION GAMA(26) ,TE 
DIMENSION FNULIM(2) ,P 
DIMENSION CR(10) ,DR(1 


EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 


ROG) sc 
(11,9) ,€ 
(ALFA(L,1 
(BETA(1,1 
(BETA(1,5 
DATA ELIM1,ELIM2, TOL 


DATA PP(1)/8.72909153 
PP(3)/1.24578576 


TOLS=LN(1.E-3) 
DATA TOLS 


DATA UPOL(1) ,CON1,CON 
6, 6666666666667E- ol, 
1.0416656566667E-01/ 


DATA RTWO,PDF,RTTP,PI 
7. 8539816339745E-01, 


DATA FNULIMN(1)/100./, 


- A NON-FATAL ERROR 


vALFA(26,4) 
“(C2 (22 
A2(52) 


(ALFA(1,3), 


(BETA(1,3) 


INDICATOR 
NORMAL RETURN - COMPUTATION COMPLETED 
X 15 LESS THAN 0.0 
ALPHA IS LESS THAN 0.0 
N 1S LESS THAN 1 

0 LAST NZ COMPONENTS OF Y SET T0 0.9 
BECAUSE OF UNDERFLOW 


GUMENTS - A FATAL ERROR 
(NZ.6T.0) 


RX ,DTM, DEN 
BETA (26,5) 


BL (52) ,B2(52) ,B3(26) 

MP (3) ,KMAX(5) ,AR(8) ,BR(10) 
P(4) 

0) 


wUPOL (10) 


1(1)) 
2(1)) 
),A1(1)) 
PAN) 
BHA) 
»B21)) 
),B3(1)) 

i Obds ; 44, f 1.E-15 
PP(2)/2.6569393226503E-01/, 


Q35SSE+00/, 
E PP(4)/7.70133574743039E-04/ 


3555 
B6559E-O1/, 


/-6.9077552799821E+00/ 


2,CONS ,CONS48 / L  OOTOOOODODODOE SO: 
SB, 355333333333 3E-01, 1.4142155623731E+00, 
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BESJ2Z11C 
BESJ211D 
BESJ211E 
BESJ211F 
BESJ2116 
BESJ211H 
BESJ2111 
BESJ211J 
BESJ 211k 
BESJ211L 
BESU211M 
BESJ211N 
BESJ2110 
BESJ211P 
BESJ2110 
BESJ211R 
BESJ2115 
BESJ2t2 

BESUZ13A 
BESJ213B 
BESJ213C 


moon 


M4 poe ~O CON O&O OWS OP re ao cn & C4 


4 Poe 


Beir O ON oO On SS 


B2(45)/ 4. 5853948516536E-05/, 82(46)/ 4.0144551389149E-05/, 
B2(47)/ 3, 5048173003133E-05/, B2(48)/ 3.0515799503435E-05/, 
B2(49)/ 2.64935611995052E-05/, B2(50)/ 2.2936363369100E-05/, 
B2(51)/ 1.9789305666402E-05/, B2(52)/ 1.7009198463541E-05/ 
DATA BS(1) / 7.3646581057258E-04/, B3(2) / 8.72790805144619E-04/, 
B3(3) / 6.22614862573514E-04/, B3(4) / 2.8599815419430E-04/, 
B3(S) / 3.8473757287937E-06/, B3(6) /-1.8790600363597E-04/, 
B3(7) /-2.9760364659456E-047, B3(8) /-3.4599812683267E-04/, 
BS3(9) /-3.53558247091604E-04/, B3(10)/-3.3571563577505E-04/, 
B3(11)/-3.04352112478904E-04/, B3(12)/-2,.6672272304761E-04/, 
BS(13)/-2.2755421412282E-04/, B3(14)/-1.8992261195456E-04/, 
B3(15)/-1.5505891859909E-04/, B3(16)/-1.2377824076187E-04/, 
B3(17)/-9.6292614771764E-05/, BS(18)/-7.2517832771442E-05/, 
B3(19)/-5. 2207002889563E-05/, #3(20)/-3.5034775051190E-G5/, 
B3(21)/-2.0848976103555E-05/, B3(22)/-8.7010609684977E-04/, 
B3(23)/ 1.1569868667510E-G6/, B3(24)/ 9.1642647412278E-06/, 
B3(25)/ 1.5647778542887E-05/, B3(26)/ 2.0822362948247E-05/ 
DATA GAMA(1L) /6.29960592494744E-01/, 
GAMA(2) /2.5198420997898E-01/, GAMA(S) /1.54790300415466E-01/, 
GAMA(4) /1.10713062414616E-01/, GAMA(S) /8.5730939552740E-02/, 
GAMA(6) /6.9716131695868E-02/, GAMA(7) /5.8608567199371E-02/, 
GAMA(8) /5.04698873535631E-02/, GAMAI9) /4.42600580689146E-02/, 
GAMA(10)/3.9372066154351E-02/, GANA(11)/3.5428319592446E-02/, 
GAMA(12)/3.2181885750210E-02/, GAMA(13)/2.9464624079116E-02/, 
GAMA (14) /2.7158167711293E-02/, GAMA(15)/2.5176827297386E-02/, 
GAMA (16) /2.3457075530608E-02/, GAMA(17)/2.1950839013491E-02/ 
GAMA (18) /2.0621082823565E-02/, GAMA(19)/1.94388240897988E-02/, 
GAMA(20)/1.838106S380068E-02/, GAMA(21)/1.7429321323196E- 02/, 
GAMA (22) /1.65685835778661E-02/, GAMA(23) /1.5796528598792E-02/, 
GAMNA(24) /1.5972950149410E- -02/, GAMA(25) /1.4419325083994E-02/, 


GAMA(26)/1. 


S8184805973534E-02/ 


TEST INPUT ARGUMENTS 


NZ= 

KT=1 

IF(N-1) 92 

Kale? 

NN=N 

IF (X) 

IF (ALPHA) 

POE 

IF (N.EQ. 1) 

Ie? 

60 TO 118 

Til 

DOMit9) l=11).N 
(IN 0F 

RETURN 

CONTINUE 

IF (ALPHA.LT.O.) 


108, 


109 


93,110,120 
91,114,116 


RETURN 


60 T0 91 
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BESJ213D 
BESJ213E 
BESJ213F 
BESJ2136 
BESJ214 
BESJ215A 
BESJ215B 
BESJ215C 
BESJ215D 
BESJ215E 
BESJ219F 
BESJ2156 
BESJ215H 
BESJ2151 
BESJ213d 
BESJ 215K 
BESJ215L 
BESJ215M 
BESJ2164 
BESJ217A 
BESJ217B 
BESJ2170C 
BESJ217D 
BESJ217E 
BESJ217F 
BESJ2176 
BESJ217H 
BESJ2171 
BESJ217J 
BESJ 217K 
BESJ217L 
BESJ217M 
BESJ217N 
BESJ218 
BESJ219 
BESU220 
BESJ221 
BESJ222 
BES ae 
BESJ2 
Rie 
BESJ22 
BESJ22 
BESJ22 
BESJ229 
BESJ23 
BESJ231 
BESJ232 
seer 
ee 
BESJ23 
BESJ237 
BESJ23 
BESJ239 


[os al oe a oe a oe 


oO 


(52) f3) (cal (se) fs) 


850 


840 


880 


DFN=DBLE (FLOAT (N) )+DBLE (ALPHA) -1.D+0 
FNU=DFN 

XO2=X*.5 

SX02=X02*X02 


DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X 
TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE 


APPLIED. 


IF(SXO02.LE. (FNUt+1.)) GO TO 850 
TA=AMAX1 (20. ,FNU) 
IF(X.GT.TA) GO TO 880 
IF(X.GT.12.) GO TO 840 
XO2L=ALOG (X02) 
NS=SX02-FNU 

60 TO 852 

FN=FNU 

FNPI=FN+1. 

XO2L=AL0G (X02) 

IS=KT 

IF(X.LE.90.5) GO TO 134 
NS=0 
DFN=DFN+DBLE (FLOAT (NS) ) 
FN=DFN 

FNPL=FN+1. 

IS=KT 

IF(N-1+NS.G6T.0) I5=3 

60 TO 134 
NS=AMAX1(36.-FNU,O.) 
DFN=DFN+DBLE (FLOAT (NS) ) 
FN=DFN 

IS=KT 

IF(N-1+NS.G6T.0) I5=3 

60 TO 130 

CONTINUE 

RTX=SQRT(X) 
TAU=RTWO#RTX 
TA=TAUtFNULIM (KT) 
IF(FNU.LE.TA) GO TO 500 
FN=FNU 

IS=KT 


UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY 


CONTINUE 

XX=X/FN 

W2=1.-XX#XX 

ABW2=ABS (W2) 

RA=SQRT (ABW2) 
IF(ABW2.67.9.2775) GO TO 200 


CASES NEAR X=FN, ABS(1.-(X/FN)¥*#2).LE.0.2775 
COEFFICIENTS OF ASYMPTOTIC EXPANSION BY SERIES 
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PESJ240 
BESJ241 
BESJ242 
BESJ243 
BESJ244 
BESJ245 
BESJ246 
BESJ247 
BESJ248 
BESJ249 
BESJ290 
BESJ251 
BESJ252 


BESJ258 
BESJ209 
BESJ260 
BESJ261 
BESJ262 
BESJ263 
BESJ264 
BESJ265 
BESJ256 
BESJ267 
BESJ268 
BESJ269 
BESJ270 
BESJ271 
BESJ272 
BESJ2753 
BESJ274 
HESJ275 
BESJ276 
BESJ277 
BESJ278 
BESJ279 
BESJ280 
BESJ2981 
BESJ282 
BESJ283 
BESJ284 
BESJ285 
BESJ286 
BESJ287 
BESJ288 
BESJ289 
BESJ270 
BESJ291 
BESJ292 
BESJ293 
BESJ294 


Iona 


ba 
bh 


ioe) noonn 


ra 


ZETA AND TRUNCATION FOR A(ZETA) AND B(ZETA) SERIES BESJ295 


BESJG296 
KMAX IS TRUNCATION INDEX FOR A(ZETA) AND B(ZETA) SERIES=NAX(2,SA) BESJ297 
BESJ298 
SA=0. BESJ299 
IF (ABW2.EQ@.G.) GO TO 21 BESJ300 
SA=TOLS/ALOG (ABW2) BESJS301 
SB=SA BESJ302 
DO 22 1=1,5 BESJ303 
KMAX (IT) =AMAX1(SA,2.) BESJ304 
SA=SA+SB BESJS05 
CONTINUE BESJ306 
KB=KMAX (9) BESJ307 
KLAST=KB-1 BESJ308 
SA=GAMA (KB) BESJ309 
DO 24 K=1,KLAST BESJ310 
KB=KB-i BESJS11 
SA=SA#W2+GAMNA (KB) BESJ312 
CONTINUE BESJ313 
Z=W2*SA BESJ314 
AZ=ABS(Z) BESJS15 
RTZ=SQRT (AZ) BESJ3516 
FNIZ=FN¥*CON2 BESJ317 
RTARY=RTZ#FN13 BESJ318 
ARY=-RTARY*#RTARY BESJ319 
AZS2=AZ#RTZ*CON1 BESJS2 
ACZ=FN*¥AZ32 BESJ321 
MEN oLE55)) (0) Wo) 227/ BESJ522 
BESJ323 
TEST FOR UNDERFLOW, 1.E-2B0=EXP(-444.), ONE WORD LENGTH BESJ324 
UP FROM UNDERFLOW LIMIT OF CDC 46600 BESJ325 
BESJ326 
IF(ACZ.GT.ELIM2) GO TO 180 BESJ327 
ARY=-ARY BESJS2 
PHI=SQRT(SQRT (SA+SA+SA+SA) ) BESJ32 
BESJ330 
B(ZETA) FOR S=0 BESJ331 
BESJ332 
KB=KMAX (3) BESJ333 
KLAST=KB-1 BESJ334 
SB=BETA(KB, 1) BESJS39 
DO 23 K=1,KLAST BESJ334 
KB=KB-1 BESJ337 
SB=SB#W2+KETA(KB, 1) BESJ358 
CONTINUE BESJ339 
KSP1=1 BESJ340 
FN2=FN#FN BESJ341 
RFEN2=1./FN2 BESJ342 
RDEN=1. BESJ343 
ASUM=1. BESJ 344 
RELB=TOL*ABS (SB) BESJ345 
BSUM=SB BESJ346 
DO 25 KS=1,4 BESJ347 
KSPI=KSP1+1 BESJ348 
RDEN=RDEN#RFN2 BESJ349 
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A(ZETA) AND B(ZETA) FOR S=1,2,35,4 


KB=KMAX (5-KS) 
KLAST=KB-1 
SA=ALFA(KB,KS) 
SB=BETA(KB,KSP1) 

DO 26 K=1,KLAST 
KB=KB-1 
SA=SA#W2+ALFA(KB,KS) 
SB=SB#W2+BETA(KB,KSP1) 
CONTINUE 

TA=SA#RDEN 

TB=SB*#RDEN 
ASUM=ASUM+TA 
BSUM=BSUM+TB 
IF(ABS(TA).LE.TOL.AND.ABS(TB).LE.RELB) GO TO 152 
CONTINUE 

CONTINUE 

BSUM=BSUM/ (FN*#FN1S) 

60 TO 400 


CONTINUE 
TAU=1./RA 
T2=1./W2 
IF(W2.GE.0.) 60 TO 30 


CASES FOR (X/FN).GT.SQRT(1.27735) 


AZ32=ABS (RA-ATAN(RA) ) 
ACZ=AZ32#FN 

CZ=-ACZ 

Z732=1,.5*AZ32 
RTZ=Z52**CON2 
FNIS=FN#*CON2 
RTARY=RTZ#FN1S 
ARY=-RTARY#RTARY 

60 T0 150 

CONTINUE 


CASES FOR (X/FN).LT.SQRT (0.7225) 
AZS2=ABS(ALOG((1.+RA)/XX) -RA) 


TEST FOR UNDERFLOW, 1.E-280 = EXP(-644.), ONE WORD LENGTH 
UP FROM UNDERFLOW LIMIT OF CDC 4600 


ACZ=AZ32*FN 

CZ=ACZ 

IF(ACZ.GT.ELIMN2) GO TO i980 
Z32=1.5#AZ32 

RTZ=Z232**CON2 
FNIS=FN**CON2 
RTARY=RTZ#FN15 
ARY=RTARY*RTARY 
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BESJ3550 
BESJ351 
BESJ352 
BESJ353 
BESJ354 
BESJ355 
BESJ356 
BESJ397 
BESJ398 
BESJ359 
BESJ340 
BESJ341 
BESJ342 
BESJ563 
BESJ364 
BESJ365 
BESUS466 
BESJ367 
BESJ548 
BESJ349 
BESJS70 
BESJ371 
BESJ372 
BESJ373 
BESJ3574 
BESJ375 
BESJ376 
BESJ377 
BESJ378 
BESJ379 
BESJ380 
BESJ381 
BESJ382 
BESJ383 
BESJ3584 
BESJ385 
BESJ386 
BESJ387 
BESJ388 
BESJS89 
BESJS90 
BESJ391 
BESJ392 
BESJ393 
BESJ394 
BESJ399 
BESJ396 
BESJ397 
BESJS98 
BESJ399 
BESJ460 
BESJ401 
BESJ402 
BESJ403 
BESJ404 


oO 


150 


101 


155 
145 


400 


402 


CONTINUE 
PHI=SQRT((RTZ+RTZ) *TAU) 
TB=1. 

ASUM=1. 

TFN=TAU/FN 
UPOL(2)=(C(1,1)*T24+C(2,1)) #TFN 
RCZ=CONI/C2Z 
CRZ32=CONS48*RCZ 
BSUM=UPOL (2) +CRZ32 
RELB=TOL¥ABS (BSUM) 
AP=TFN 

KS=0 

KP1=2 

RZDEN=RCZ 

DO 155 LR=2,8,2 


COMPUTE TWO U POLYNOMIALS FOR NEXT A(ZETA) AND B(ZETA) 


LRPI=LR+1 
DO 101 K=LR,LRP1 

KS=KS+1 

KPL=KP1+1 

S1=C(1,K) 

DO 102 J=2,KP1 
S1=S1*T2+C(J,K) 

CONTINUE 

AP=AP*TFN 

UPOL(KP1)=AP*S1 

CR(KS) =BR(KS) #RZDEN 
RZDEN=RZDEN*RCZ 
DR(KS)=AR(KS) #RZDEN 
CONTINUE 

SUMA=UPOL (LRP1) 
SUMB=UPOL (LR+2) +UPOL (LRP1) #CRZ32 
JU=LRP1 

DO 151 JR=1,LR 

JU=JU-1 
SUMA=SUMA+CR (JR) #UPOL (JU) 
SUMB=SUNB+DR(JR) #UPOL (JU) 
CONTINUE 

TB=-TB 

IF(W2.GT.0.) TB=ABS(TB) 
ASUM=ASUM+SUMA*TB 
BSUM=BSUM+SUMB#TB 

IF (ABS(SUMA).LE. TOL. AND. ABS(SUMNB).LE.RELB) GO TO 165 
CONTINUE 

TB=RTARY 

IF(W2.6T.0.) TB=-TB 
BSUM=BSUN/TH 


CONTINUE 
CALL JAIRY(ARY,RTARY,ACZ,AI,DAI) 
TEMP (1S)=PHI* (AI *ASUM+DAI*BSUM) /FNIS 
GO TO (401,202,650), IS 

TEMP (1) =TEMP (3) 
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BESJ405 
BESJ406 
BESJ407 
BESJ408 
BESJ409 
BESJ410 
BESJ411 
BESJ412 
BESJ413 
BESJ414 
BESJ4195 
BESJ416 
BESJ417 
BESJ418 
BESJ419 
BESJ 420 
BESJ421 
BESJ422 
BESJ423 
BESJ424 
BESJ425 
BESJ426 
BESJ427 
BESJ428 
BESJ429 
BESJ430 
BESJ431 
BESJ432 
BESJ433 
BESJ434 
BESJ4355 
BESJ434 
BESJ437 
BESJ438 
BESJ439 
BESJ440 
BESJ441 
BESJ442 
BESJ443 
BESJ444 
BESJ445 
BESJ446 
BESJ447 
BESJ448 
BESJ449 
BESJ450 
BESJ451 
BESJ452 
BESJ4535 
BESJ454 
BESJ455 
BESJ456 
BESJ4397 
BESJ458 
BESJ499 


oO 


oa 


ET=1 BESJ440 


401 I8=2 BESJ441 
DFN=DFN-1.D+0 BESJ462 
FN=DEN BESJ443 
60 TO 130 BESI444 

BESJ445 
SERIES FOR (X/2)#**2.LE.NU+1 BESJ446 
BESJ447 

134 CONTINUE BESJ448 
GLN=GAMLN (FNP1) BESJ449 
ARG=FN*XO2L-GLN BESJ470 
IF(ARG.LT.-ELIM1) GO TO 123 BESJ471 
EARG=EXP (ARG) BESJ472 

300 CONTINUE BESJ473 
§sil. BESJ474 
AK=3. BESJ475 
eh. BESJ474 
Tain BESJ477 
S1=FN BESJ478 
MO 123 Reilly BESJ479 
$2=T2+51 BESJ480 
T=-T#5X02/52 BESJ481 
S=S+T BESJ482 
IF(ABS(T).LT.TOL) GO TO 127 BESJ495 
T2=T2+AK BESJ484 
AK=AK+2. BESJ485 
S1=S1+FN BESJ486 

125 CONTINUE BESJ487 

7 CONTINUE BESJ488 
TEMP (1S) =S*EARG BESJ499 
60 TO (301,202,600), IS BESJ490 

301 EARG=EARG*FN/ X02 BESJ491 
DFN=DFN-1.D+0 BESJ492 
FN=DFN BESJ493 
1S=2 BESJ494 
GO TO 360 BESJ495 

BESJ496 
SET UNDERFLOW VALUE AND UPDATE PARAMETERS BESJ497 
BESJ498 

180 Y(NN)=0, BESJ499 
NN=NN-1 BESJ500 
DFN=DFN-1.D+0 BESJSO1 
FN=DFN BESJ502 
IF (NN-1) 170,171,130 BESJ 503 

17 KT=2 BESJ504 
15=2 BESJ505 
GO TO 130 BESJ504 

123 Y(NN)=0. BESJS07 
NN=NN-1 BESJSO8 
FNPL=FN BESJ509 
DFN=DFN-1.D+0 BESJS10 
FN=DFN BESJ511 
IF(NN-1) 170,172,173 BESJS12 

IW2MKeT=2 BESJ5S13 
1S=2 BESJS14 
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300 


ott 
512 


IF(SX02.LE.FNP1) GO TO i33 
60 TO 130 
ARG=ARG-XO2L+ALOG (FNPI1) 
IF(ARG.LT.-ELIM1) GO TO 123 
60 TO 134 

NZ=N-NN 

RETURN 


BACKWARD RECURSION SECTION 


CONTINUE 

NZ=N-NN 

IF(KT.EQ@.2) GO TO 250 
CONTINUE 

BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA 
¥Y<(NN) =TEMP (1) 

Y(NN-1) =TEMP (2) 
IF(NN.EQ@.2) RETURN 
DX=X 

TRX=2.D+0/DX 
DTIM=DFN#TRX 

TM=DTM 

K=NN+1 

DO 230 I[=3,NN 

K=K-1 
¥(K=-2) =TM#Y (K-21) -Y¥ CK) 
DTM=DTM-TRX 

TM=DTM 

CONTINUE 

RETURN 

Y(1)=TEMP (2) 

RETURN 


ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN 
OSCILLATORY REGION X.GT.MAX(29, NU), PROVIDED THE LAST MEMBER 


OF THE SEQUENCE IS ALSO IN THE REGION. 


CONTINUE 
IN=ALPHA-TAU+2. 
IF(IN.LE.O) GO TO 562 


INPI=IN+1 
DALPHA=ALPHA-FLOAT(INP1) 
KT=1 

GO TO 5i1 

DALPHA=ALPHA 

IN=06 

IS=KT 


ARG=X-PIDT#DALPHA-PDF 
SA=SIN(ARG) 

SB=COS (ARG) 
RA=RTTP/RTX 

ETX=8. #X 

DX=DALPHA 

DX=DX+DX 

DTM=DX*DX 
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BESJ568 
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ao 
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T2=DTM-1.D+0 
T2=T2/ETX 
$2=T2 
RELB=TOL*ABS(T2) 
TL=ETX 

S1=1. 

FN=1, 

AK=8. 

DO 504 K=1,13 
TI=TI+ETX 
FN=FN+AK 
DX=FN 
TRX=DTM-DX 
AP=TRX 
T2=-T2*AP/T1 
S1=S1+T2 
TLSTI+ETX 
AK=AK+8. 
FN=FN+AK 
DX=FN 
TRX=DTH-DX 
AP=TRX 

T2= T2*AP/T1 
§2=52+T2 
IF(ABS(T2).LE.RELB) GO TO 505 
AK=AK+8. 


4 CONTINUE 
3 TEMP(1S)=RA*(S1*SB-S2*SA) 


60 TO (506,507) ,1S 
DALPHA=DALPHA+1. 
15=2 

TB=SA 

SA=-SB 

SB=TB 

60 TO 303 


FORWARD RECURSION SECTION 


IF(KT.E8.2) GO TO 250 
S1=TEMP (1) 

S2=TEMP (2) 

TX=2./X 

TM=DALPHA*TX 
IF(IN.E@.0) GO TO 520 


FORWARD RECUR TO INDEX ALPHA 


DO S10 I=1,1N 

=52 

S2=TM*52-S1 

TM=TM+TX 

S$1=5 

CONTINUE 

IF(NN.EQ@.i1) GO TO 535 


$=82 
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BESJ604 
BESJ605 
BESJ4606 
BESJ407 
BESJG608 
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BESJ410 
BESJ411 
BESJé12 
BESJé13 
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920 


600 


6950 


677 


603 


S2=TH#S2-S1 
TM=TH+TX 
51=5 
CONTINUE 


FORWARD RECUR FROM INDEX ALPHA TO ALPHA+N-1 


Y(t) =S1 

¥(2)=52 

IF(NN.E@.2) RETURN 
DO 530 I=3,NN 
Y¥(T)=TM*¥Y (1-1) -¥ (1-2) 
TM=TM+TX 

CONTINUE 

RETURN 

Y¥(1)=S2 

RETURN 


BACKWARD RECURSION WITH NORMALIZATION BY 
ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES. 


CONTINUE 

COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION 
KM=AMAX1(S.-FN,O.) 

TFN=FN+FLOAT (KM) 


TA=X02L-TA 

TB=-(1.-1.5/TFN)/TFN 
IN=CE/(-TA+SQRT (TAX TA-CE*¥TH))+1.5 
IN=IN+KN 

GO TO 603 

CONTINUE 


COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION 


GLN=AZ32+RA 

IF(ARY.GT.30.) GO TO 475 
RDEN=(PF(4) #ARY+PP(3)) *ARY+1. 
RZDEN=PP (1) +PP (2) #ARY 
TA=RZDEN/RDEN 

IF(W2.LT.G.16) GO TO 451 
TB=GLN/RTARY 

60 TG 677 

TB=(1. 259992104940. 16798947 30#W2) /FNLS 
GO 10 677 

CONTINUE 

TA=CONL*¥TCE/ACZ 

TA=((0, O4935827160#TA-O. LLLLLLLL11) #TAtG. 6666666667) €TAXARY 
IF(W2.LT.0.10) GO TO 451 
TB=GLN/RTARY 

IN=TA/TB+1.5 

IFC IN. GT.INLIM) GO TO 4062 
DX=FLOATC(IN) 

DTM=DFN+DX 

DX=X 

TRX=2.D+0/Dk 

DTM=DTM#TRX 
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BESJ627 
BESJ628 
BESJ429 
BESJ630 
BESJ631 
BESJ652 
BESJ433 
BESJ634 
BESJ6355 
BESJ436 
BESJ637 
BESJ438 
BESJ639 
BESJ4640 
BESJ641 
BESJ642 
BESJ643 
BESJ644 
BESJ649 
BESJ4646 
BESJ647 
BESJ648 
BESJ649 
BESJ4650 
BESJ651 
BESJ652 
BESJ695 
BESJ654 
BESJ655 
BESJ696 
BESJ4&57 
RESJ658 
BESJ659 
BESJ460 
BESJ661 
BESJ642 
BESJ6635 
BESJ664 
BESJ6465 
BESJ466 
BESJ667 
BESJ468 
BESJ469 
BESJ670 
BESJ471 
BESJ6/72 
BESJ473 
BESJ4674 
BESJ4735 
BESJ676 
BESJ677 
BESJ478 
BESJ479 


aod 


a 


605 


601 


604 
615 


oy 


TM=DTM 
TA=0. 
TB=TOL 
KK=1 
CONTINUE 


BACKWARD RECUR UNINDEXED 


DO 601 I=1,1N 

S=TH 

TB=TM*TE-TA 

TA=S 

DTM=DTM-TRX 

TM=DTM 

CONTINUE 
NORMALIZATION 
IF(KE.NE.1) GO TO 4604 
TA=(TA/TB) #TEMP (3) 
TB=TEMP (3) 

KK=2 

IN=NS 

IF(NS.NE.0) GO TO 605 
Y (NN) =TB 

NZ=N-NN 
IF(NN.E@.1) RETURN 
5=TB 

TB=TM#TB-TA 

TA=S 

DTM=DTM-TRX 

TM=DTM 

K=NN-1 

Y(K)=TB 

TF (NN.EQ@.2) RETURN 
KM=K-1 


BACKWARD RECUR INDEXED 


DO 602 I=1,KM 
Y(K-1) =TM#Y (K) -¥ (K#1) 
DTM=DTM-TRX 

TM=DTM 

K=K-1 

CONTINUE 

RETURN 


CONTINUE 
NZ=-2 
RETURN 
CONTINUE 
NZ=-3 
RETURN 
CONTINUE 
NZ=-1 
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ioe) 


CDC 6600 ROUTINE 


AND ITS DERIVATIVE DAI(X) FOR JBESS 
INPUT: X - ARGUMENT, COMPUTED BY JBESS, X UNRESTRICTED 
RX — RX=SORT(ABS(X)), COMPUTED BY JBESS 
C - C=2.*#(ABS(X)*#1.5)/3., COMPUTED BY JBESS 
OUTPUT: AI - VALUE OF FUNCTION AI(X) 
DAI - VALUE OF THE DERIVATIVE DAI(X) 
WRITTEN BY D.E. AMOS, S.L. DANIEL & M.K WESTON 
DIMENSION AKI (14) ,AK2(23) ,AKS (14) 
DIMENSION AJP(19) sAUN(19) ,AC15) ,B(15) 
DIMENSION DAK1 (14) ,DAK2(24) ,DAK3 (14) 
DIMENSION DAJP(19) ,DAIN(19) ,DA(15) ,DB(15) 
DATA Ni,N2,N5,N4/14,23,19,15/ 
DATA M1,M2,M3,M4/12,21, 17,13 / 
DATA FP112,CON1,CON2,CONS,CON4 ,CONS/ 


1 


5 
< 


Cl roe 
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RETURN 
END 


SUBROUTINE JAIRY(X,RX,C,Al,DAI) 


1.5089969389958E+00, 
3.8000458986729E-01, 


DATA AKI(1) 
AK1(2) /-1. 
AK1(4) / 8. 
AK1(6) / 1. 
AK1(B) /-1 


DATA AK2(1) 
AK2(2) / 
AK2(4) / 
AK2(6) / 
AK2(8) / 
AK2(10)/ 
AK2(12)/ 
AK2(14)/ 
AK2(16)/ 
AK2(18)/ 
AK2(20)/ 
AK2(22)7 

DATA AK3(1) 
AKS (2) 
AKS(4) 
AK3 (4) 
AKS (8) 


AK3(14)/ 
DATA AJP(1) 
AJP (2) 


98762 
»8729020974102E- 06/, 
»21135041695591E-07/, 

29245470528 265E- 097, 
»26863505659574E-10/, 
~JOL41035901371E-11/, 
©§4464450473454E-12/, 

»7010603055935E- 
-1977479916481E-14/, 

144734079E- 


et et oO ee lO) Rl ed oo CO I | 


HON 
fe 
cal ike 
Laie 
AK3(10)/-6. 
AR S(L2)/-1. 


fHIkq 


RASC] 


JAIRY COMPUTES THE AIRY FUNCTION AI(X) 


/ 2. 2042309098779E- 
2929024278770E-01/7, 
2284415200554E-04/, 
6382428017212E-05/7, 


»29621999359335E-07/, 
AKL(190)/ 1. 
AKL(12)/ 2. 
AK1(14)/-35. 


Rese see ae VASES g 
O32535257423563E-12/, 
16169897 78508E- Maid 
/ 2.7436615086960E- 


3979096975 F690E-03/, 


2924875E- o4/, 
1713 1890E-05/, 


13/, 


/ 2.8027 
7812704284438E-03/, 
b3524996526900E-06/, 
D2 29435502291 6E-09/, 
2440825180026E-11/, 
S63752444639351E-135/, 
1270513469106E-14/, 


-2.46480324312435E-16/ 
4 7.7895295643758E- 


B43565654568GE-01/, 
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01/, 

AKL (3) 
AK1 (3) 
AK1(7) 
AK1 (9) 
AKI (11)/- 


SS 


1. 
20 
3. 
8. 
3. 


AKI(13)/-1. 


O1/, 
AK2 (3) 
AK2(5) 
AK2(7) 
AK 2 (9) 


AK2(13)/ 
AK2(15)/-1. 
AK2(17)/-9, 


AK2(19)/-4. 


AK2(21)/-4, 
AK2(23)/-3. 
OLA 
AK3(3) / 
AK3(5) / 
ARSC) if 
AK3(9) / 
AK3(11)/ 
AK3(13)/ 


02/, 
AIP(S) / 


fio 
(No 
JI 5 
7-4. 
AK2(11)/-3. 
-2,0390244716791E- 


— cO cn cn ~O 


9.03515471619678E+00, 
8. 6602540378444E-01/ 


OSBBLL6335919E-02/, 
34614345891 23E-04/ 

0690258957319E-07/, 
2290815882567E-09/, 
3916546561568E-11/, 
1067954609788E-14/, 


5733922062119E-03/, 
1212491739992E-04/, 
3680422537055E-06/, 
7589279396229E-07/, 
0924537427061E-08/, 
09/~" 
3670476763957E-10/, 
3139829654843E- 12/, 
4384026199096E-13/, 
5076010450328E-14/, 
1907704086507E-15/ 


© OS542257962900E-05/, 
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»4713840457655E-10/7, 
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AIRY39 
AITRY40 
AIRY41 
AIRY42 
AIRY43 
AITRY44 
AIRY45 
AIRY46 
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AIRYS3 


oOo ON Oo OF & hd by 
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AJP(4) / 3.0534272427761E-02/, AJP(S) /-4.9542470251308E-03/, 
AJP (4) /-1.7274955256395E-03/, AIP(7) / 2.4313763783919E-04/, 
AJP(8) / 5.0456477751708E-05/, AJP(9) /-6.16316582469521E-06/, 
AJP (10) /-9.0398574551077E-07/, AIP(11)/ 9.7024377835588E-08/, 
AJP(12)/ 1.0963945330520E-08/, AJP(13)/-1.0471633058877E-09/, 
AJP (14) /-9.6035944134465E-11/, AJP(15)/ §.2535878945413E-12/, 
AJP (16)/ 6.3612343901877E-13/, AIP(17)/-4.9662961411602E- 14/, 
AJP (18) /-3.2981028892962E-15/, AJP(19)/ 2.3579825203110E-16/ 
DATA AJN(1) / 3.8049788761724E-02/, 
AJN(2) /-2.4531954184555E-01/, AJN(3) / 1.6582062370270E-01/, 
AIN(4) / 7.4933004581879E-02/, AJN(S) /-2.6347628810664E-02/, 
AIN(5) /-5.9253559730498E-03/, AIN(7) / 1.4474440958980E-03/, 
AIN(8) / 2.1831183132222 -04/, AIN(9) /-4.1065207748030E-05/, 
AIN(10) /-4.6687499417177E-06/, AIN(L1)/ 7.1521880727714E-07/, 
AIN(12)/ 6.5296477085463E-08/, AIN(13)/-8.4428402756595E-09/, 
AJN(14) /-6.4419615897698E-10/, AIN(15)/ 7.20802 28650528E-11/ 
AIN(14)/ 4.7246543171785E-12/, AJN(17) /-4.6602263254704E- 137, 
AIN(18) /-2.6776271038919E-14/, AJN(19)/ 2, 3616131657002E-15 
DATA A(1) / 4.9027542474279E- D1/, A(2) / 1 .5764727794620E-03/, 
A(3) /-9.6619596314031E-05/, A(4) / 1.3591608028882E-07/, 
A(5) / 2.9815734265486E-07/, Ald) /-1.868247675599B8E-08/, 
A(7) /-1.0368573766714E-09/, A(B) / 3.28660818434335E-10/, 
A(9) /-2.570914106327BE-11/, A(10)/-2.32357655300468E-12/, 
A(11)/ 9.5752327904826E-13/, A(12)/-1.2034082864972E-13/, 
A(13) /-2.9090771677072E-15/, A(14)/ 4.5565445458015E-15/, 
A(15) /-9.9900387481026E-16/ 
DATA B(1) / 2.7859355280308E-01/, Bi2) /-3. Seiteey ieee 53E-03 
B(3) /-2.3114967738499E-05/, B(4) / 4.7131784226356E- nee 
B(S) /-1.1241590793133E-07/, B(6) /-2. NOLOOSDL AREAS, 
B(7) / 2.6094807530219E-09/, B(8) /-3.5509813610122E-11/, 
B(9) /-3.508499784238B8E-11/, B(10)/ 5.8300718795420E-12/, 
B(11)/-2.0464482875333E-13/, H(12)/-1.1052917947674E-13/, 
B(13)/ 2.9772477803878E-14/, B(14)/-2.8820511100994E-15/, 
B(15)/-3.3265631169617E-16/ 
DATA NID,N2D,N3D,N4D/14,24,19,15/ 
DATA M1D,M2D,M3D,M4D/12,22,17,13/ 
DATA DAK1(1) / 2.0456784230789E-01/, 
DAKi(2) /-6.6132273990566E-02/, DAKI(3) /-8.4984590098929E-03/, 
DAKL(4) / 3.1218349155629E-03/, DAKL(S) /-2.7001648982943E-04/, 
DAK1(5) /-6.3563629867939E- 08/4 DAKI(7) / 3.0239771240951E-06/, 
DAK1(8) /-2.1831119533009E-07/, DAKI(9) /-5.3619428933293E-10/, 
DAKIi(10)/ 1.1309803562231E- ae DAK1 (11) /-7.4302383462907E-11/, 
DAK1(12)/ 4.2890417082689E-13/, DAKL{13)/ 2.2381092575454E-13/, 
DAKI (14) /-1, 3914013564118E- 14/ 
DATA DAK2(1) / 2.9333234388323E-01/, 
DAK2(2) /-8.0619878474311E-03 3/, DAK2(3) / 2.42540172333514E-03/, 
DAK2(4) /-.8229754885024E-04/, DAK2(5) / ines 8642775118E-04/, 
DAK2(6) /-4.9745744768406E-05/, DAK2(7) / 1.3209068123950E-05/, 
DAK2(8) /-3.4952824044494E-06/, DAK2(9) / 9.2436245107894E-07/, 
DAK2(10)/-2.4473257152187E-07/, DAK2(11)/ 6.4930783764891E-08/, 
DAK 2 (12) /-1.7271762150154E-09/, DAK2(13)/ 4.6072576360466E-09/, 
DAK2 (14) /-1.2324905529155E-09/, DAK2(15)/ 3.3062040949810E-10/, 
DAK2(15)/-8.8925209977240E-11/, DAK2(17)/ 2.3977351987930E-11/, 
DAK2 (19) /-S.4801392115345E-12/, DAK2(19)/ 1.7551013202373E-12/, 
DAK2(20) /-4. 7630382983344E-13/, DAK2(21)/ i. 29498241 10081E-13/, 
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DAK2 (22) /-3.5267 
-2.62786914342 
/ 2.8467582881135E-01/7, 


DAK2(24)/ 
DATA DAKS(1 
DAKS(2) / 
DAKS(4) 7 
DAKS(6) / 
DAKS(8) / 
DAKS(10)/ 
DAK3S (12)/ 
DAKS(14)/ 
DATA DAJP(1 
DAJP(2) /- 
DAJP(4) / 
DAJP(6) /- 
DAJP(8) / 
DAJP (10) /- 
DAJP(12)/ 


DAJP(14)/-3., 


DAJP(16)/ 
DAJP(18)/ 
DATA DAUN(1 
DAJN(2) / 
DAJN(4) /- 
DAIN(4) / 
DAJN(8) /- 
DAJN(1G)/ 
DAJN(12)/- 
DAJN(14)/ 
DAJN(16) /- 
DAIN(18) / 
DATA DAI) 

DA(3) 


221043E-14/, 
29E-15/ 


»9507307261908E-03/, 
»84907283944634E-06/, 
»0592563445715E-69/, 
~ 563557688851 354E-11/, 
»8857435378444E-13/, 
»1737476261721E-14/, 
»0957477309706E-16/ 
) ff 
1.20262933468882E-01/, 
1.6794842923050E-02/, 
8.4556029509887E-04/, 
2.2982786094548E-05/, 
3.76343991134692E-07/, 
4, 29611335 200301E-09/7, 
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DAJP (3S) 
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6200315158592E-15/, 


8348115033798E-05/, 
0141849117858E-07/, 
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39631765351 04E-12/, 
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11164930827496-03/, 
617697761721 4E- 
9729580465652E-08/, 
4475282664204E-10/, 
9965506584722E-i2/, 
6209895256874E-13/, 
5938481128449E-16/, 
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0 
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4421 2835341992E-05/, 
SBO40S51841871E-04/, 
2428689470978E-09/, 
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1867422109358E-12/, 
9OL9057660871E-14/, 
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DO 105 1=1,h1 
J=J-1 

TEMPL=F 1 
FI=TT¥F1-F2+AK1 (J) 
F2=TEMP1 

CONTINUE 
AL=T#F1-F2+AK1 (1) 


J=N1D 

F1=DAK1 (J) 
F2=0. 

DO 106 I=1,M1D 
J=J-1 

TEMP1=F1 


FL=TT¥F1-F2+DAK1 (J) 


F2=TEMP1 
CONTINUE 


DAI=-(T#F1-F2+DAK1 (1) ) 


RETURN 


CONTINUE 
T=(X+X-CON2) ¥CONS 
We Sul kee L 

J=N2 

FL=AK2(J) 

F2=0. 

pO 155 I=1,M2 
J=d-1 

TEMP1=Fl 
FL=TT#FI-F2+AK2 (3) 
F2=TEMP 1 

CONTINUE 
RTRX=SORT (RX) 
EC=EXP(-C) 


AI=EC¥(T#FI-F2+AK2 


J=N2D 
FL=DAK2(d) 
F2=0. 

DO 156 1=1,M2D 
J=J-1 

TEMP1=F 1 


(1))/RTRX 


FI=TT¥FI-F24DAK2 (0) 


F2=TEMP 1 
CONTINUE 


DAL=-EC#(T#F1-F2+DAK2(1)) #RTRX 


RETURN 


CONTINUE 
T=10. /C-1. 
TT=T+T 

J=N1 
FI=AK3(J) 
F2=0. 

DO 205 I=1,M1 
J=J-1 
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AIRY164 
AIRY165 
AIRY166 
AIRY147 
AIRY168 
AIRY149 
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AIRY171 
AIRY172 
AIRY173 
AIRY174 
AIRY175 
AIRY176 
AIRY177 
AIRY178 
AIRY179 
AIRY180 
AIRY181 
AIRY182 
AIRY183 
AIRY184 
AIRY185 
AIRY186 
AIRY187 
AIRY188 
AIRY189 
AIRY190 
AIRY191 
AIRY192 
AIRY193S 
AIRY194 
AIRY195 
AIRY196 
AIRY197 
AIRYi98 
AIRY199 
AIRY200 
AIRY201 
AIRY202 
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AIRY204 
AIRY205 
ALRY2Z06 
AIRY 207 
AIRY208 
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AIRY213 
AIRY214 
AIRY215 
AIRY216 
AIRY217 
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ba 
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cn 


300 


309 


306 


TEMP1=F 1 

FL=TT#F1-F2+AK3 (J) 

F2=TEMP1 

CONTINUE 

RTRX=SQRT (RX) 

EC=EXP(-C) 
AI=EC#(THFI-F2+AKS(1))/RTRK 
J=N1D 

FL=DAKS(J) 

F2=0, 

DO 206 I=1,MiD 

J=J-1 

TEMP1=F1 
FL=TT#F1-F2+DAKS (J) 
F2=TEMPI 

CONTINUE 
DAI=-RTRX#EC#(T#FI-F2+DAKS(1)) 
RETURN 


CONTINUE 
IF(C.G67T.5.) 60 TO 350 
T=.44#C-1. 

TT=T+T 

J=N3 

FL=AJP (J) 
EL=AUN(J) 

F2=0. 

E2=0. 

DO 305 I=1,M5 
J=d-1 

TEMP1=F 1 

TEMP2=E1 
FLSTT#Fi-F2+AUP (J) 
EL=TT#E1-E2+AUN (J) 
F2=TEMPI 

E2=TEMP2 

CONTINUE 


AL=(T#EL-E2+AIN(L))-X¥(T#FL-F2+ANP (1) ) 


J=N3D 

FL=DAUP (J) 
EL=DAIN(J) 

F2=0. 

E2=0. 

DO 306 I=1,M3D 

J=J-1 

TEMP1=F 1 

TEMP2=E1 

Fl = TT¥FL-F2+DAJP (J) 
El= TT#E1-E2+DAIN(J) 
F2=TEMP1 

E2=TEMP2 

CONTINUE 


DAI=X#X*¥(T#FI-F2+DAJP(1))+(T#E1-E2+DAIN(1)) 


RETURN 
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350 


CONTINUE 
T=10./C-1. 

TT=T+T 

J=N4 

FL=A(d) 

E1=B(J) 

F2=0, 

E2=0. 

DO 310 I=1,H4 
J=d-1 

TEMP1=F 1 

TEMP2=E1 
FL=TT#F1-F2+A(d) 
E1=TT¥E1-E2+B(J) 
F2=TEMP1 

E2=TEMP2 

CONTINUE 
TEMP1=T#F1-F2+A (1) 
TEMP2=T¥EL-E2+B (1) 
RTRX=SORT (RX) 
CV=C-FPI12 
CCV=COS(CV) 
SCV=SIN(CV) 
Al=(TEMP1*¥CCV-TEMP2*SCV) /RTRX 


J=N4D 

FLi=DA(J) 
EL=DB(J) 

F2=0, 

E2=0. 

DO 3ii 1=1,M4D 
J=J-1 

TEMP1=F 1 
TEMP2=E1 


FI=TT#F1-F2+DA(J) 
EL=TT#E1-E2+DB(d) 
F2=TEMF1 

E2=TEMP2 

CONTINUE 
TEMPL=T#F1-F2+DA(1) 
TEMP2=T#E1-E2+DB (i) 
EL=COV*CONS+. S¥SCV 
E2=SCV*CONS-.S*CCV 
DAI=(TEMPL*E1-TEMP2*E2) *RTRX 
RETURN 

END 

FUNCTION GAMLN(X) 


WRITTEN BY D. E. AMOS, SEPTEMBER, 1977. 
REFERENCES 
* SAND-77-1518 


* COMPUTER APPROXIMATIONS BY J.F.HART, ET.AL., SIAM SERIES IN 
APPLIED MATHEMATICS, WILEY, 19468, P.135-134. 

* NBS HANDBOOK OF MATHEMATICAL FUNCTIONS, AMS S55, BY 
M. ABRAMOWITZ AND I.A. STEGUN, DECEMBER. 1955, P.257 
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AIRY274 
ATRY2735 
AIRY276 
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ALRY27 
AIRY279 
ATRY280 
AIRY281 
AIRY282 
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AIRY284 
AIRY285 
AIRY286 
ATRY287 
AIRY288 
AIRY289 
AIRY290 
ATRY291 
AITRY292 
AIRY293 
AIRY294 
AIRY295 
AIRY296 
AIRY297 
AITRY298 
ATRY299 
AIRY300 
AIRY3O01 
AIRY302 
AIRY303 
AIRY304 
ATRY305 
AIRY306 
AIRY307 
AIRY308 
AITRY309 
ATRY210 
ATRYS11 
AIRYS12 
ATRY313 
AIRY314 
AIRYS15 
AIRYS14 
AIRYS17 
AIRY318 
GLNI 
GLN2 
GLNS 
GLN4 
GLNS 
GLN6 
GLN7 
GLNS 
GLN9 
GLN1LO 


NIOOIAINIAAMAAWAAAABAAHANHM 


GAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR 


X.GT.G. A RATIONAL CHEBYSHEY APPROXIMATION IS USED ON 
B.LT.X.LT.1000., THE ASYMPTOTIC EXPANSION FOR X.SE.1000. 
A RATIONAL CHEBYSHEV APPROXIMATION ON 2.LT.X.LT.3. FOR 
O.LT.X.LT.8. AND X NON-INTEGRAL, FORWARD OR BACKWARD 
RECURSION FILLS IN THE INTERVALS 0.LT.X.LT.2 AND 
S.LT.X.LT.8. FOR X=1.,2 3100., GAMLN IS SET TO 
NATURAL LOGS OF FACTORIALS. 
DESCRIPTION OF ARGUMENTS 
INPUT 
X SeXGh 0 
OUTPUT 
GAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT xX 
DIMENSION GLN(100) ,P(S) ,Q(2) ,PCOE(9) ,Q@COE(4) 
DATA XLIM1,XLIM2,RTWPIL/ EG 9 aoe » 97.189385322047E- 
DATA P(1)/7.663451880000E-04/, P(2)/-5.940956105200E-04/, 
CSS) / To PESTS NOSE TSER O47, PCA) /-2.7777777546577E-03/, 


i 

2 
DATA 
DATA 


PW he 


ted Poo 


oon oon Pine Oona one 


DATA 


Cn ee tel Pa oe 


RCo ies 
Q(1)/- 
PCOE(1)/2. 


PCOE(3)/1. 
PCOE(S)/2, 
PCOE(7)/2. 
PCOE(9) /6. 


OCOE(1)/1. 
QCOE(3)/9. 


GLN(1) /0. 
GLN(4) /1., 
GLN(4) /4, 
GLN(8) /8.52 
GLN(10)/1. 
GENCI2) 7107 
GLN(14)/2.2 
GLN(16)/2, 
GLN(18)/3.3 
GLN(20)/3., 
2)/4,.5 
GLN(24)/5. 
GLN(26)/5. 
WAG 4 
GLN(30)/7. 
GLN(32)/7., 
GLN(34)/8. 
GLN(34)/9, 
Nao Fe 
GLN(40)/1., 


GLN(2 


GLN(28 


SLN(38 


BLN(42)/1. 


SLN(44)/1. 


GLN(45)/1. 


GLN(48)/1., 
"445657439463E+02/, 
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BLN(S0) /1 


edrdend 


ITIIIIIIT 


-02/ 


7777777777786 - O3/, 


973786644810E-03/ 
O931159546710E-01/ 
159943128451E+00/ 
0782472535179E+01/ 
2003583800713E+01/ 


DODDDDD00NDOE+00/, 


Q(2)/8. 3333535355333 3E- 
2581943959028E- 
980671310204E-01/7, 
338067999387E+00/, 
503677233002E+01/, 


»PCOE(2)/9. 
,PCOE(4)/3. 
»PCOE(4) /6. 
»PCOE(8)/3. 


QCOE(2 


02/ 
03/, 


)/-8.906016659498E+00/, 


822521104714E+00/ ,QCOE (4) 76. 20038380071 3E+O1/ 


O/7, GLN(2) /0.0/, 
TILTS94G9228E+00/, 
TE7AIIT4A27E2E+O0/, 
es ere 
OLB2748008E+01/, 
O2307B84587E+01/, 
216385512E+01/, 
9927138384E+01/, 
SSOS0734S014E+01/, 
FS359BB41B8720E+01/, 
38013889848E+01/, 
160667554776E+01/, 
BOG3IE0S22298E+01/, 
JI7SSB62701E+01/, 
12570289671 7E+01/, 
B09222555332E+01/, 
JOS44S701758E+01/, 
213617560369E+01/, 
P33051245479E+017, 
066317502506E+02/ 
14G342117815E+02/, 
219330815 1S4E+02/, 
29 ee eae 


3480272243 73E+ 


ee ee ies Ie 
coc Gy 
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GLN (3) 
GLN(3) 


GLN(7) /6. 
GLN(9) /1. 
GLN(I1)/1. 
GLN(13)/1. 
GLN(15)/2. 
GLN(17)/3., 


GLN(19)/ 


3.63590445 


GLN(21) 7/4, 
GLN(23)/4., 


GLN (25) 


GLN(27)/6. 
GLN(29)/6, 
GLN(31)/7. 


GLN(33)/8. 
GLN(35)/8. 
GLN(37)/9, 


GLN{39)/1. 
GLN(41)/1, 


GLN(43)/1. 
GLN(45)/1. 


GLN(47)/1. 
GLN(49)/1. 
GLNCS1)/1. 


A27 


75.931471805599E-01/, 
/3.178053830348E+00/, 


S79251212010E+00/, 
050460290274E+01/, 
S10441257308E+01/, 
9987214495466EtO1/, 
S19122118274E+01/, 
DS71BS01060BE+91/, 
Z20B03E+01/, 
255561646075E+01/, 
B471181351B84E+0i/, 


/S.A7TB4A7T2Z9S9BLIE+O1/, 


261701761 00E+01/, 
7TEB974AZLSTLISE+Oi/, 
Fees Oe e eS Eun. 

LSS7ISPASSLZE+O1/, 
898082754220E+01/, 
S71G69454214Et+01/, 
O29681996145E+02/ 
LOSZ06397148E+02/, 
177718813997E+92/, 
205172711 494E+02/, 
S295Z5750356E+02/, 
40673923648 2E+02/, 
48477766951 8E+02/, 


AND 


o1/ 


GLNi1 
GLN12 
GLNIS 
GLN14 
GLNIS 
GLN14 
GLN17 
GLN18 
GLNI9 
GLN20 
GLN21 
GLN22 
GLN23 
GLN24 
GLN25 
GLN26 
GLN27 
5LN28 
GLN29 
GLNSO 
GLN3S1 
GLNS2 
GLN33 
GLN34 
GLN3S 
GLNS6 
GLN37 
GLN38 
GLN39 
GLN40 
GLN41 
GLN42 
GLN43 
GLN44 
GLN45 
GLN46 
GLN47 
SLN48 
GLN49 
GLNSO 
GLNGi 
GLNS2 
BLNSS 
GLNS4 
6GLNSS 
GLNS6 
SLNS7 
GLNS8 
GLNO? 
GLNSO 
GLN61 
GLNé2 
GLN6S 
GLN64 
GLNGS 


a 


Lap) 


(val (sa) (o2) 


an 


10 


15 


— ~0 CO “I cr 


L4 poe OO coONM co Cn & 4 Po or 


oon or cn 


DATA 


DATA 


GLN(S2)/1. 
GLN(S4)/1. 
BLN(S4) /1 
GLN(S8)/1. 
GLN(60)/1, 
GLN(61)/1. 
GLN(S6S)/1. 
GLN(65)/2. 
GLN(67)/2. 
GLN(69) 7/2. 
GLNI71)/2. 
GLN(73)/2. 
GLN(75)/2. 
GLN(77)/2. 
GLN(79)/2. 
GLN(81)/2. 
GLN(83)/2. 
GLN(85)/2. 
GLN(S7)/3., 
GLN(89)/3. 
GLN(91)/3. 
GLN(93)/3, 
GLN(95)/3. 
GLN(97)/3. 
GLN(99)/3. 


S24095925945E+02/, 
6035511 2B215SEt+02/, 


~68S274454484E+62/, 


7639584840 70E+02/, 
8453358288614 E+G2/ 
8862917542357E+02/, 
968661814572Z9E+02/, 
O51681994826E+02/, 
LS5322414945E+02/, 
219564418191 E+92/, 
SO439043546S58E+02/, 
SB97B3895618E+02/, 
4§73729140962E+02/, 
562211355500E+02/, 
649216497985E+02/, 
736735124285 7E+02/, 
B24742926876Et+02/, 
FL S2Z39500943E+027, 
OO2Z2074B6470E+02/, 
U91L641935802E+02/, 
18152463596202E+02/, 
271852877038E+02/, 
So2511819792E+02/, 
4337940706235E+02/, 
J45359085S5194E+02/, 


NDX=X 
T=X-F 
POT 


LOAT(NDX) 
£@.6.0) 


DX=XLIMI-X% 


IF (DX. 


RATIONAL CHEBYSHEY APPROXIMATION ON 2. 


LT.0.0) 


NXM=NDX-2 
PX=PCOE (1) 

DO 10 1=2,9 
PX=T#PX+PCOE (1) 
QX=OCDE(1) 

DO 15 I=2,4 
QX=T#OX+QCOE (1) 


DGAM= 


IF (NX 
IF (NX 


BACKWARD RECURSION 


PX/QX 
M.GT.9) 


MEG. 0) 


60 TO 2 
60 T0 2 


60 TO 51 


GO TO 40 


DGAM=DGAN/(1.+T) 


IF (NX 


DGAM= 


M.EQ.-1) 
DGAN/T 


60 T 


Ac 
aud 


GAMLN=ALOG (DGARM) 
RETURN 


FORWARD RECURSION FOR 3.LT.X.LT.9 
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GLN(62)/1. 


GLN(54)/2 

GLN(646)/2. 
GLN(48)/2. 
GLN(70)/2, 


GLN(72)/2. 


GLN(74) 7/2. 
GLN(76)/2. 
GLN(78)/2. 
GLN(S0)/2. 


GLN(82)/2. 
GLN(84)/2. 
GLN(84)/2, 


SLN(88)/3. 
GLN(90)/3. 
GLN(92)/3. 
GLN(94)/3. 


GLN(96)/3. 
GLN(98)/3. 
GLN(100)/3. 


Cag) fin 
GLN(S55)/1. 
GLN(S7)/1. 
GLN(S9)/1. 
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92739047 287B8E+02/, 
OLOO9S163993E+02/, 
0934258675 25E+02/, 
L77Sh9S41140E+02/, 
2619054983235 7E+02/, 
BATOL7234428E+02/, 
AS26884 900S0E+92/ : 
S18904022097E+02/, 
SOS649F09719E+027, 
9291097 651L0E+02/, 
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26534991 267E+02/, 
5171788719695 +02 bf 
40B8150588708E+02/, 
499541180408E+02/ 
5913542053695E+02/ 


3 FOR GAMMA(X) 
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on] 


lao) 


oO 


40 


af 


XX=2.4T 

DO 24 1=1,NkR 
DGAMN=DGAN*X Xx 
XX=KX#1. 
GANLN=ALOGS(DGAN) 
RETURN 


X.GT.xXLIM1 


RX=1./% 

RXX=RX#RX 

IF((X-XLIM2).LT.0.) GO TO 41 
PX=O(1) #RXX+O(2) 
GAMLN=PX#RX+(X-.5) #ALOG(X) -K+RTWPIL 
RETURN 


X.LT.XLIM2 


PX=P (1) 
SUM=(X-.5) #ALOG (X) -X 
DO 42 1=2,5 
PX=PX#RXX+P (1) 
CONTINUE 
GAMLN=PX#RX+SUM+RTWEIL 
RETURN 


TABLE LOOK UP FOR INTEGER ARGUMENTS LESS THAN OR EQUAL 100. 


IF(NDX.GT.100) GO TO 40 
GAMLN=GLN(NDX) 

RETURN 

END 
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APPENDIX B: NOTATION 


Incident wave amplitude 
Softy 

Gravitational acceleration 

Water depth 

v=1 

Bessel function of the first kind 

Wave number 

Non-negative integers 

Polar coordinate 

Temporal coordinate 

Horizontal coordinate 

Horizontal coordinate 

Bessel function of the second kind 
Vertical coordinate 

Incident wave angle 

Free surface displacement 

Polar coordinate 

Angle related to wedge angle 

Value related to wedge angle 

Velocity potential function 

Horizontal component of the velocity potential function 
Defined in Equation 12 

Incident wave velocity potential function 
Reflected wave velocity potential function 


Scattered wave velocity potential function 


Finite cosine transform of ¢ 


Wave radian frequency 
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